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ABSTRACT 


The temperature evolution of icosahedral medium-range order formed by interpenetrating icosahedra in CuZr metallic glass¬ 
forming liquids was investigated via molecular dynamics simulations. Scaling analysis based on percolation theory was 
employed, and It Is found that the size distribution of clusters formed by the central atoms of Icosahedra at various temperatures 
follows a very good scaling law with the cluster number density scaled by and the cluster size S scaled by 11 - 
respectively. Here Tc Is scaling crossover-temperature, t and a are scaling exponents. The critical scaling behaviour suggests 
that there would be a structural phase transition manifested by percolation of locally favoured structures underlying the glass 
transition. If the liquid could be cooled slowly enough but without crystallization Intervening. Furthermore, It is revealed that 
when icosahedral short-range order (ISRO) extends to medium-range length scale by connection, the atomic configurations of 
ISROs will be optimized from distorted ones towards more regular ones gradually, which significantly lowers the energies of 
ISROs and introduces geometric frustration simultaneously. Both factors make key Impacts on the drastic dynamic slowdown of 
supercooled liquids. Our findings provide direct structure-property relationship for understanding the nature of glass transition. 


Introduction 

The nature of glass and glass transition is the deepest and most interesting unsolved problem in the solid state science'^. With 
a generic definition, numerous systems spanning a broad range of length scale such as atomic and colloidal systems, foams, and 
granular materials, can be considered as glass when certain conditions are satisfied^. Metallic glass, as a relatively “simple” 
glassy system for scientific research of glass transition and promising industrial materiaP, has attracted much attention and 
interest of scientists from broad research fields®’^. 

The CuZr metallic glass-former has been extensively investigated*"'®. It exhibits good glass-forming ability, and its binary 
chemical composition reduces the complexity of local atomic structures, which makes this system a good model for the study of 
structure-property relationships in liquids and glasses. Previous studies demonstrate that icosahedral short-range order (ISRO) 
is closely correlated with the slow dynamics and dynamical heterogeneity during glass formation*"'^’ It has also been 
demonstrated that the formation of icosahedral medium-range structures via the interpenetration of ISROs plays a key role in the 
dynamic slowdown'"’ ", The more the ISROs are connected, the slower the dynamics of the connected ISROs is. This implies 
that the icosahedral medium-range order formed by the percolation of ISROs may be related to the glass transition. So far, the 
concept of percolation of ISROs has been widely used for understanding the relationship between stmctural evolution and glass 
transition in metallic glass-forming liquids. However, the question is that no specific percolation theory or scaling analysis 
has been derived to establish a direct link between the percolation of ISROs and glass transition in the metallic glass-forming 
liquids. In addition, the population of ISROs in CuZr metallic glass-forming liquids sometimes is not very high so that ISROs 
do not even percolate as glass transition occurs'^. For example, CusoZrso is a good glass-former in both experiments'* and 
computer simulations'*’ '*. However, the fraction of ISROs in CusoZrso at 300 K is found to be less than 4% in simulation". 
Even in the inherent structure of CusoZrso, it is less than I5%®. To establish the link between percolation of ISROs and glass 
transition, the nearest-neighbor atoms of ISROs were also taken into account for the percolation. Meanwhile, the percolation 
of ISROs in CuZr metallic glass-forming liquids could be influenced by system size in simulation, leading to uncertainty for 
understanding the strcumre-dynamics relationships in these systems. Therefore, it is highly desirable to develop a scaling 
analysis based on percolation theory to establish a quantitative description of percolation of ISROs and link to glass transition. 
Moreover, although ISROs are found to correlate with the slow dynamics and the connection of ISROs makes it even slower, 
the physical origin is still not very clear. 


In this work, molecular dynamics simulations were performed for CusoZrso with realistic interatomic potential (see 
Methods). We investigated the connection of ISROs via interpenetrating or volume-sharing. Graph theory was introduced 
to characterize the the clusters formed by the central atoms of interpenetrating ISROs at different temperamres as the CuZr 
metallic glass-forming liquids are cooled down, and the equilibrium cluster size distribution was analyzed. Scaling analysis 
based on percolation theory^^’'^^ was conducted for the cluster size distribution. It is found that the cluster size distributions at 
various temperamres collapse together and follow a good scaling law, as the cluster size is S scaled by 11 — T’c/T’p^^'^ and the 
cluster number density is scaled by 5^^, respectively. The scaling analysis suggests that there could exist a geometric phase 
transition of percolation of locally favoured stmctures, once the metallic glass-forming liquids are quenched slowly enough (in 
the hypothetical limit of infinitely long relaxation time) but without crystallization intervening, and glass transition may be 
related to the percolation of locally favoured structures. Furthermore, it is revealed that as ISROs are connected together, the 
atomic conhgurations of connected ISROs are optimized towards more regular icosahedra. The optimized ISROs enhance the 
geometric frustration, and the local energies are signihcantly lowered, which stabilizes the structure of liquids and slows down 
the dynamics in glass transition. 

Results and Discussions 

Scaling analysis for the percolation of icosahedral network 
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Figure 1 . Number density distributions of clusters formed by the connection of ISROs via volume sharing at different 
temperatures during cooling. The inset illustrates a cluster with size S = 2 formed by two volume-shared ISROs. 

First, we analyzed the size distributions of clusters formed by the connection of central atoms of the icosahedra at various 
temperatures. In our analysis, graph theory was applied for the construction of connection of icosahedra and the resulting 
icosahedral network (see Methods). As shown in Fig. 1, the size distributions decrease monotonically as cluster size increases 
and exhibit similar behaviour at different temperatures. The size distributions in small size part follow power-law behaviour, but 
deviate in larger size part, which signals there is a hnite characteristic cluster size in the system and, for a given temperature T, 
the characteristic cluster size marks the crossover between a power-law behaviour and a rapid decay of n{S,T), qualitatively. As 
temperature decreases, more and more larger clusters are formed in the supercooled liquids. It would be expected that the cluster 
size distribution could asymptotically approach a pure power-law behaviour as temperature further decreases and approaches 
some critical point if possible. The similarity of the clusters size distributions at different temperatures is reminiscent of a 
general scaling behaviour. Therefore, to explore the scaling law behind the cluster size distributions, a general scaling ansatz*® 
for the cluster number density n{S,T) based on percolation theory was employed. That is, 

n{S,T)<^S-^fiS/Sc), (1) 


( 2 ) 



where the characteristic cluster size Sc diverges as a power-law in terms of the distance of T from Tc'. 
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Figure 2. The scaled size distribution follows an excellent scaling behaviour with the scaled cluster size St- Here S is the 
cluster size. Tc is the critical point. T and a are scaling critical exponents. It can be seen that all the cluster number densities 
fall onto the same curve representing the graph of scaling function f{S/Sc)- 


and the function / is known as the scaling function for the cluster number density. Generally, the expression of / varies from 
system to system, and dimension to dimension. Analytic solution of / is non-trivial in most cases. However, the asymptotic 
behaviour of / can provide sufficient information of the scaling behaviour of the size distributions. According to percolation 
theory, for 5/5c <C 1, the scaling function is approximately constant, and for S/Sc S> 1 it decays rapidly. To our knowledge^®, 
except for one-dimension percolation, this behaviour of / is quite universal in the scaling analysis based on above scaling 
ansatz. 

To reveal the underlying scaling behaviour of the cluster size distributions, following above scaling ansatz, we scaled 
the cluster number density n{S,T) with S^'^ and the cluster size S with |1 — respectively. Figure 2 shows that 

the scaled cluster size distributions at different temperatures collapse onto a master curve representing the graph of scaling 
function f{S/Sc)- Here Tc=700 K is the scaling crossover-temperature. t = 2.05 and (1=0.4 are the scaling exponents. These two 
exponents are different from the values obtained in three-dimensional site percolation^because the ISROs percolate in a 
continuous space, but not on a discrete periodic lattice. As shown in Fig. 2, the scaling function f{S/Sc) decays rapidly as 
the scaled cluster size Sj ^ I, which indicates that there is a characteristic cluster size in the system. The evolution of the 
characteristic cluster size Sc shows a divergence behaviour as temperature is approaching the critical point Tc- AtT = Tc, the 
characteristic cluster becomes infinite, so that n{S,T) S^'^ f{S/Sc) becomes 

n(S,T)o.S-\ (3) 

because it can be seen from Fig. 2 that the scaling function f{S/Sc) will approach a non-zero constant for St I- The 
scale-free behaviour of the cluster size distribution at threshold shows some information about the geometric properties of 
the percolating cluster, indicating that the incipient infinite cluster has an internal fractal geometry. The fractal dimension 
df ^ 2.86 of the structures formed by connected ISROs at the critical point can be calculated from the scaling relation of 

T = d/df + l, (4) 

where (/ = 3 is the spatial dimension*^. It can be seen that the value of df is not that small, which indicates that the atomic 
packing of incipient infinite cluster is not so loose. On the other hand, just as the characteristic cluster size Sc diverges as T is 
approaching Tc, the correlation length associated with the connected ISROs also diverges. For a particular characteristic cluster 
size Sc, the associated radius of gyration defines a characteristic length scale that is proportional to the correlation length 
According to percolation theory, one has 

(5) 

so that the temperature dependence (as T Tc) of the correlation length ^ can be characterized as 

^oc\1-Tc/T\-\ ( 6 ) 
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which shows a power-law behaviour with a divergence at T = Tc. The corresponding critical exponent v = 0.875 can be 
determined by the scaling relation^^ of 

V = l/cJii/'= (t- l)/(J(i. (7) 

As shown above, in percolation, once the cluster number density is known, all other quantities can be derived. All results 
indicate that the percolation of ISROs forming the so-called icosahedral medium-range orders indeed correlates with the glass 
transition in CuZr metallic glass-forming liquids. We argue that df of the incipient infinite cluster emerging at Tc cannot be 
calculated directly by the box-counting method. The reason is as follows: Since the crossover temperature Tc is blow the glass 
transition temperature (over 900 K) obtained in the previous study^^, the atomic structures at 700 K obtained in simulations 
are in non-equilibrium states and sensitive to the cooling-rate. As a result, dj obtained from box-counting analysis for the 
atomic structures at 700 K will be cooling-rate dependent. This is not exactly the same df derived in percolation theory for the 
equilibrium structural phase transition. Therefore, df obtained from the universal scaling relation (t = d/df + 1) is generic. 
It is also worth noting that Tc is not cooling-rate dependent, because it is a critical temperature derived from the equilibrium 
structural phase transition of the modelling system, and all data for determining Tc are generated from the relatively equilibrium 
supercooled liquid states. This is similar to To in Vogel-Fulcher-Tamman (VFT)^’^ equation. 

Physical origin of dynamic slowdown associated with percolation 
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Figure 3. The self-intermediate scattering functions of the central atoms of ISROs with different node degree k in CuZr 
metallic glass-forming liquid at 1000 K. The inset shows an exponential dependence of scaled relaxation times Xkj'^a on k/12. 

The above scaling analysis demonstrates that the percolation of ISROs in CuZr metallic glass former during cooling is 
closely related to glass transition. It is still not clear why the percolation of ISROs could contribute to slowing down of 
dynamics. As mentioned above, the connectivity of ISROs signihcantly influences the dynamical property of local structures'*^. 
In order to emphasize the importance of the medium-range structures formed by the connection of ISROs to the dynamic 
slowdown, we calculated the self-intermediate scattering functions^^ (SISFs) of the central atoms of icosahedra with different 
node degree 

1 / A'r \ 

( Eexp{/^- [c(f) -r;(0)]} \ , (8) 

where the sum is over all central atoms with degree k, rjit) is the location of atom j at time t. q is chosen with the amplitude 
approximately equal to the value of the first peak position in the static structure factor, and (•) denotes the ensemble average. 
The relaxation time Xk is determined by F^{q, Xk) = Figure 3 shows the SISFs as a function of time for central atoms 
of the icosahedra with different k values (due to the very rare connection of k > 5, SISFs of k > 5 were not calculated***.) in 
supercooled metallic glass-forming liquid at 1000 K. It can be seen that all SISFs for different k values exhibit non-exponential 
decay behaviour, and the decay becomes much slower as k value increases, dramatically depending on k. The inset in Fig. 3 
explicitly shows that the scaled relaxation time Xk/Xa increases exponentially with the value of k/12 {Xa is the average structural 
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relaxation time, and 12 is the maximum value of node degree)^**. It can be seen that the relaxation time of the atoms with large 
k is even more than 10 times of the Ta- This finding reveals that the medium-range structures formed by the connection of 
ISROs fundamentally influences the relaxation dynamics of local atomic structures. As the icosahedral medium-range order 
percolates during quenching, the resulting dynamic slowdown spreads in the whole system, leading to the sluggish dynamics, 
which Anally contributes to glass transition. 

It is also very interesting to investigate how the icosahedral medium-range order influences the atomic symmetry of 
icosahedra themselves. All icosahedra in the MD modeled samples are actually distorted from the ideal icosahedron, because 
from a geometrical viewpoint, icosahedral clusters cannot All the entire three-dimensional space without partially breaking of 
the five-fold rotational symmetry^®. As observed by Frank over half a century ago, the ideal icosahedral arrangement indeed 
has a significantly lower energy than fee atomic arrangement^^. The question is whether the connection of ISROs will promote 
dense packing and five-fold local symmetry of ISROs. If so, the self-aggregation effect of icosahedra^ will naturally tend to 
minimize the local energy density, slow the atomic dynamics, and lead to great geometric frustration. In order to get deep insight 
into the above discussion, we analyzed the local atomic symmetry of ISROs and its dependence on the connectivity degree k. 
To analyze the local atomic symmetry of ISROs, the bond orientational order (BOO) parameter introduced by Steinhardt et 
was adopted, in which the BOO of the £-fold symmetry is defined as a 2f -f 1 vector; 


where is spherical harmonics and A,- is the number of bonds of atom i with its nearest neighbor atoms. In the analysis, one 
use the rotational invariants defined as 


qi = 
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where the term in brackets in the invariants (11) is the Wigner 3 j symbol^®. For fee, bcc, and hep symmetry, the value of is 
close to zero, while for perfect icosahedron, W(, = —0.169757, so that is sensitive to reflect the five-fold atomic symmetry 
features of ISROs in metallic glass-forming liquids. Figure 4 shows that the average values of IFg of ISROs are not -0.169757, 
but in between -0.06 and -0.10, suggesting that the ISROs in metallic glass-forming liquids are not perfect, but distorted with 
partially fee symmetry, in agreement with the previous studies^®. Furthermore, as shown in Fig. 4, with increasing node degree 
k, W(, decreases to more negative values. This indicates that as the node degree increases, on the range we studied, the atomic 
configurations of ISROs are optimized towards more regular icosahedra. The configuration optimization of ISROs produces 
significant geometric frustration in the supercooled liquids. 

We also calculated the formation energy Ep^c of ISROs in the supercooled metallic liquids and investigated the dependence 
of the formation energy on node degree k. The formation energy of ISROs can be calculated by the formula of 

^ ~ ^ref,a) , ( 13 ) 


where Ep j is the potential energy of the y'th atom in an icosahedron and is the reference energy of the element 
In our calculations, the chemical potentials of Cu and Zr in crystal structures (fee for Cu and hep for Zr) were used as the 
reference energies for Cu an Zr atoms, respectively. Therefore, in equation (13), the effect of the chemical compositions in an 
icosahedron was eliminated, and the formation energies of ISROs with different chemical compositions can be comparable. 
Thus, we can analyze the formation energies of ISROs with different degree k. The inset histogram in Fig. 4 clearly shows 
that with increasing degree k, the formation energy of ISROs decreases monotonically, indicating that the formation of the 
icosahedral medium-range order leads to lower energies, and the structures become more stable. All of above results indicate 
that the connection between different icosahedra will indeed generate a positive feedback for optimization of the local ISROs. 
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Figure 4. The dependence of W(, of an ISRO on its node degree k in supercooled liquids at 900 K (black) and 1000 K (red), 
respectively. The inset histogram shows the change of the formation energies of ISROs with different node degree k. Two 
atomic configurations of ISROs with k = 0 and k = A obtained from MD simulation were also presented for geometric 
comparison. 


This effect will minimize the local energy density, slow down the atomic dynamics, cause great geometric frustration, and 
finally contributes to glass transition. 

In summary, we carried out MD simulations for CuZr metallic glass-forming liquids to investigate the relationship between 
percolation of ISROs and glass transition. As the system is supercooled, the ISROs tend to connect with each other to form big 
clusters. The cluster size distributions at different temperatures follows an excellent scaling law, which indicates that the cluster 
size distribution evolves toward a power-law behaviour, and an infinite size cluster is formed as the system is approaching a 
critical point. Therefore, there would be a structural phase transition manifested by percolation of locally favoured structures 
underlying the glass transition if the liquid could be cooled slowly enough but without crystallization intervening. Furthermore, 
the percolation of ISROs and the formation of medium-range orders optimize the atomic configurations of ISROs towards 
more regular icosahedra, enhancing the geometric frustration and minimizing the local energy density, which leads to the 
dynamic slowdown of the metallic glass-forming liquids. Our findings suggest that the geometric phase transition manifested 
by percolation of locally favoured structures could be critical for the understanding of the nature of glass transitions. 


Methods 

Molecular dynamic simulations 

In this work, molecular dynamic (MD) simulations were carried out for CusoZrso metallic alloy using the LAMMPS package^ ^ 
The interatomic interaction was described by the well-developed embedded-atom method potential for CuZr alloys^^. The 
sample contains 10000 atoms being randomly distributed in a cubic box with periodic boundary condition applied in three 
dimensions. The MD step is 2 fs. At first, the sample was equilibrated at 2000 K for 4 ns (2,000,000 MD steps) in NPT (P=0) 
ensemble with Nose-Hoover thermostat and barostat. The liquid was then quenched at a rate of 1 K/ps down to its target 
temperature (1600K, 1500 K, 1400 K, 1300 K, 1200 K, 1100 K, and 1000 K, respectively). During cooling, the box size 
was adjusted to maintain zero pressure. At each temperature, the atomic configuration was relaxed in NPT (P=0) ensemble 
for another 2 ns (1,000,000 MD steps) for the analysis of physical properties (500 atomic configurations were collected for 
ensemble average). The structural relaxation time (Ta) of our modelling system at 1000 K was order of magnitude of 10 ps, so 
the ensemble average window for analysis was much longer than Xa of the system at each temperature we studied. Note that at 
temperatures higher than the glass transition point, the structures corresponding to (metastable) equilibrium can be achieved 
quite rapidly. Therefore, in general, the effect of cooling rate on the analysis of structural, dynamic, or other physical properties 
of our modelling system is eliminated. 
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Icosahedral network and node degree 

The local atomic structures in supercooled liquid samples at different temperatures were analyzed by the Voronoi tessellation 
method^^"^^ and identified in terms of the Voronoi index {n^,n 4 , 115 , 11 ( 1 ), where n,(/=3,4,5,6) denotes the number of i-edged 
faces of a Voronoi polyhedron. In our analysis, a cutoff distance of 5 A was chosen so that the Voronoi index distribution was 
converged. To characterize the connectivity of ISROs in the system, we introduced the graph theory^^~^^. In our scheme***’^®, 
the central atom (with Voronoi index (0,0,12,0)) of an icosahedron is treated as a node, and two nodes are considered to be 
connected if they are the nearest neighbors with each other, that is, the two icosahedra are interpenetrating or volume-sharing. 
The choice of connection criterion is reasonable based on recent experimental observation^**, in which the extent of icosahedral 
short-range order to form medium-range order is consistent with a facing-sharing or interpenetrating configurations. Therefore, 
the case of interpenetrating configuration was considered in our analysis. Property variation depending on the connection 
criterion will be discussed in detail in future works. With the definitions of nodes and edges abstracted from the atomic 
modelling system, we established the icosahedral network. Based on the scenario of graph theory, the node degree k was 
defined as the number of other nodes directly connected to it. The maximum value of node degree for the central atom of an 
icosahedron is 12. Our results suggest that node degree introduced from graph theory is a good nonlocal (to some extent) order 
parameter for classifying atoms with distinct properties. 
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